home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / colnew.sci < prev    next >
Text File  |  1999-09-16  |  5KB  |  175 lines
  1. function [z,z1]=col1()
  2. fixpnt=0
  3. m=[4]
  4. ncomp=1
  5. aleft=1
  6. aright=2
  7. zeta=[1,1,2,2]
  8. ipar=0*ones(1,11)
  9. ipar(3)=1
  10. ipar(4)=2
  11. ipar(5)=2000
  12. ipar(6)=200
  13. ipar(7)=1
  14. ltol=[1,3]
  15. tol=[1.e-11,1.e-11]
  16. res=aleft:0.1:aright
  17. z=colnew(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
  18.     fsub,dfsub,gsub,dgsub,guess)
  19. z1=[]
  20. for x=res,z1=[z1,trusol(x)]; end;
  21.  
  22. function [df]=dfsub(x,z)
  23. df=[0,0,-6/x**2,-6/x]
  24.  
  25. function [f]=fsub(x,z)
  26. f=(1 -6*x**2*z(4)-6*x*z(3))/x**3
  27.  
  28. function [g]=gsub(i,z)
  29. g=[z(1),z(3),z(1),z(3)]
  30. g=g(i)
  31.  
  32. function [dg]=dgsub(i,z)
  33. dg=[1,0,0,0;0,0,1,0;1,0,0,0;0,0,1,0]
  34. dg=dg(i,:)
  35.  
  36. function [z,mpar]=guess(x)
  37. // unused here
  38. z=0
  39. mpar=0
  40.  
  41. function [u]=trusol(x)
  42.   u=0*ones(4,1)
  43.       u(1) = .25* (10.*log(2.)-3.) * (1.-x) +0.5* (1./x+ (3.+x)*log(x) - x)
  44.       u(2) = -.25* (10.*log(2.) - 3.) + .5 * (-1./x/x + log(x) + (3.+x)/x - 1.)
  45.       u(3) = .5 * (2./x**3 + 1./x -3./x/x)
  46.       u(4) = .5 * (-6./x**4 - 1./x/x + 6./x**3)
  47.  
  48.  
  49. function [z,zf]=col2(flag)
  50. ipar=0*ones(1,11)
  51. //     define constants, print a heading.
  52.       gamma = 1.1d0
  53.       eps = .01d0
  54.       dmu = eps
  55.       eps4mu = eps**4/dmu
  56.       xt = sqrt(2.d0*(gamma-1.d0)/gamma)
  57. //   define no. of differential equations.
  58. fixpnt=0
  59. iflag=0
  60. ncomp = 2
  61. //   orders
  62. m=[2,2]
  63. //   interval ends
  64.       aleft = 0.d0
  65.       aright = 1.d0
  66. //   locations of side conditions
  67. zeta=0*ones(1,4)
  68.       zeta(1) = 0.d0
  69.       zeta(2) = 0.d0
  70.       zeta(3) = 1.d0
  71.       zeta(4) = 1.d0
  72. //   ipar  values
  73. //   a nonlinear problem
  74.       ipar(1) = 1
  75. //   4 collocation points per subinterval
  76.       ipar(2) = 4
  77. //   initial uniform mesh of 10 subintervals
  78.       ipar(3) = 10
  79.       ipar(8) = 0
  80. //   dimension of real work array  fspace  is 40000
  81.       ipar(5) = 40000
  82. //   dimension of integer work array  ispace  is 2500
  83.       ipar(6) = 2500
  84. //   (these dimensions of  fspace  and  ispace
  85. //    enable  colsys  to use meshes of up to 192 intervals.)
  86. //   print full output.
  87.       ipar(7) = 1
  88. //   initial approximation for nonlinear iteration is provided
  89. //   in  solutn
  90.       ipar(9) = 1
  91. //   a regular problem
  92.       ipar(10) = 0
  93. //   no fixed points in the mesh
  94.       ipar(11) = 0
  95. //   tolerances on  all components
  96.       ipar(4) = 4
  97.       ltol=1:4
  98.       tol=1.e-5*ones(1,4)
  99. //   call  colsys
  100. res=aleft:0.05:aright
  101. if flag=1 then 
  102. z=colnew(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
  103.     fsub1,dfsub1,gsub1,dgsub1,guess1)
  104. else 
  105. z=colnew(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
  106.     'f','df','g','dg','gu')
  107. end
  108. zf=[   0.00        0.00000D+00    0.20414D+01    0.10225D-31   -0.90397D+00
  109.        0.05        0.10207D+00    0.20414D+01   -0.45265D-01   -0.90794D+00
  110.        0.10        0.20414D+00    0.20414D+01   -0.90926D-01   -0.91982D+00
  111.        0.15        0.30621D+00    0.20414D+01   -0.13738D+00   -0.93962D+00
  112.        0.20        0.40828D+00    0.20414D+01   -0.18502D+00   -0.96735D+00
  113.        0.25        0.51035D+00    0.20415D+01   -0.23425D+00   -0.10030D+01
  114.        0.30        0.61245D+00    0.20430D+01   -0.28545D+00   -0.10466D+01
  115.        0.35        0.71528D+00    0.21066D+01   -0.33905D+00   -0.10992D+01
  116.        0.40        0.83213D+00   -0.24518D+00   -0.39612D+00   -0.12054D+01
  117.        0.45        0.17751D-01   -0.35755D+01   -0.44540D+00   -0.71354D+00
  118.        0.50        0.22512D-01    0.12361D+00   -0.48036D+00   -0.67007D+00
  119.        0.55        0.25869D-01    0.48526D-01   -0.51169D+00   -0.58108D+00
  120.        0.60        0.28099D-01    0.41111D-01   -0.53831D+00   -0.48234D+00
  121.        0.65        0.29911D-01    0.29912D-01   -0.55981D+00   -0.37641D+00
  122.        0.70        0.30874D-01    0.55520D-02   -0.57587D+00   -0.26595D+00
  123.        0.75        0.30032D-01   -0.45165D-01   -0.58642D+00   -0.15667D+00
  124.        0.80        0.25524D-01   -0.14662D+00   -0.59175D+00   -0.60454D-01
  125.        0.85        0.13751D-01   -0.34695D+00   -0.59307D+00   -0.14009D-02
  126.        0.90       -0.12516D-01   -0.75283D+00   -0.59330D+00   -0.28623D-01
  127.        0.95       -0.69427D-01   -0.16508D+01   -0.59906D+00   -0.24811D+00
  128.        1.00        0.26514D-13    0.11926D+03   -0.62542D+00   -0.88763D+00 ]
  129. z=z'
  130. zf=zf(:,2:5)
  131.  
  132. function [z,dmval]=guess1(x)
  133.       cons = gamma * x * (1.d0-.5d0*x*x)
  134.       dcons = gamma * (1.d0 - 1.5d0*x*x)
  135.       d2cons = -3.d0 * gamma * x
  136.       if x > xt then   z=[0,0,-cons,-dcons]
  137.       dmval(2) = -d2cons ;
  138.       else z=[ 2.d0 * x, 2 , -2*x+cons,-2+dcons]
  139.       dmval(2) = d2cons ;
  140.       end
  141.       dmval(1) = 0.d0
  142.  
  143.  
  144. function [f]=fsub1 (x,z)
  145.       f=[0,0]
  146.       f(1) = z(1)/x/x - z(2)/x + (z(1) - z(3)*(1.d0-z(1)/x) - ...
  147.     gamma*x*(1.d0-x*x/2.)) / eps4mu
  148.       f(2) = z(3)/x/x - z(4)/x + z(1)*(1.d0-z(1)/2.d0/x) / dmu
  149.  
  150.  
  151. function [df]=dfsub1 (x, z)
  152.       df=0*ones(2,4)
  153.       df(1,1) = 1.d0/x/x +(1.d0 + z(3)/x) / eps4mu
  154.       df(1,2) = -1.d0/x
  155.       df(1,3) = -(1.d0-z(1)/x) / eps4mu
  156.       df(1,4) = 0.d0
  157.       df(2,1) = (1.d0 - z(1)/x) / dmu
  158.       df(2,2) = 0.d0
  159.       df(2,3) = 1.d0/x/x
  160.       df(2,4) = -1.d0/x
  161.  
  162. function [g]=gsub1(i, z)
  163.      g=[z(1),z(3),z(1), z(4)-0.3d0*z(3)+0.7d0]
  164.      g=g(i)
  165.  
  166. function [dg]=dgsub1 (i, z)
  167.      dg=0*ones(4,1)
  168.      select i 
  169.     case 1 then dg(1)=1
  170.     case 2 then dg(3)=1
  171.     case 3 then dg(1)=1
  172.     case 4 then dg(4)=1,dg(3)=-0.3
  173.      end
  174.  
  175.